%% get the GDP
clc; clear all;
[val, text, raw] = xlsread('../data/VNIBOR.xlsx','GDP');
%% transform data to make variance constant (using log or box-cox transfrom)

gdp = val(:,1);
% lgdp = log(gdp);
[lgdp, lambdas] = boxcox(gdp);
dates = datenum(text(2:end,1));

figure(4)
subplot(2,1,1); plot(dates,gdp); title('GDP'); datetick('x',12);
subplot(2,1,2); plot(dates,lgdp); title('log GDP'); datetick('x',12);

%% 
figure(2)
subplot(2,1,1);plot(lgdp);
lgdp4s = lgdp(5:end)-lgdp(1:end-4);
subplot(2,1,2);plot(lgdp4s);
title('Log GDP');

figure(3)
subplot(2,2,1);autocorr(lgdp);
subplot(2,2,2);parcorr(lgdp);

subplot(2,2,3);autocorr(lgdp4s);
subplot(2,2,4);parcorr(lgdp4s);
hold off
%%

model = arima('Constant',0,'D',1,'Seasonality',4,...
              'MALags',1,'SMALags',4);
          
          
% Y0 = Y(1:5);
[fit,VarCov] = estimate(model,lgdp);


M = 20;
N = length(lgdp);
[Yf,YMSE] = forecast(fit,M,'Y0',lgdp);
exp(Yf(1:2))
upper = Yf + 1.96*sqrt(YMSE);
lower = Yf - 1.96*sqrt(YMSE);

figure(6)
plot(lgdp,'Color',[.75,.75,.75])
hold on
h1 = plot(N+1:N+M,Yf,'r','LineWidth',2);
h2 = plot(N+1:N+M,upper,'k--','LineWidth',1.5);
plot(N+1:N+M,lower,'k--','LineWidth',1.5)
xlim([0,N+M])
title('Forecast and 95% Forecast Interval')
legend([h1,h2],'Forecast','95% Interval','Location','NorthWest')
hold off